#include <iostream>
#include <fstream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLine.h"

using namespace std;

void CalcMuonFlux(){

  //float Ratio = 163/7128.;

  ifstream fin("./MuonFlux.txt");
  if(!fin){std::cout<<"!!! No Input File !!!"<<std::endl;return;};

  TH1F *hSpectrumMu_p=new TH1F("hSpectrumMu_p","hSpectrumMu_p",1e6,1.0,1.0e5);
  TH1F *hSpectrumMu_n=new TH1F("hSpectrumMu_n","hSpectrumMu_n",1e6,1.0,1.0e5);
  TH1F *hSpectrumMu  =new TH1F("hSpectrumMu"  ,"hSpectrumMu"  ,1e6,1.0,1.0e5);

  //--------------- Read Data
  int counter=0;
  while( !fin.eof() ){
    float Energy,neutronF,protonF,alphaF,mu_pF,mu_nF,electronF,positronF,photon;
    fin >> Energy >> neutronF >> protonF >> alphaF >> mu_pF >> mu_nF >> electronF >> positronF>> photon;
    hSpectrumMu_p->Fill(Energy,mu_pF);
    hSpectrumMu_n->Fill(Energy,mu_nF);
    counter++;
  }

  //--------------- Set All bin
  int nbin=hSpectrumMu_p->GetNbinsX();
  for(int ibin=0;ibin<nbin;ibin++){
    float vbin_p=hSpectrumMu_p->GetBinContent(ibin+1);
    float vbin_n=hSpectrumMu_n->GetBinContent(ibin+1);
    if( ibin>0 ){
      if( vbin_p==0. )	hSpectrumMu_p->SetBinContent(ibin+1,hSpectrumMu_p->GetBinContent(ibin));
      if( vbin_n==0. )	hSpectrumMu_n->SetBinContent(ibin+1,hSpectrumMu_n->GetBinContent(ibin));
    }
    hSpectrumMu->SetBinContent(ibin,hSpectrumMu_p->GetBinContent(ibin)+hSpectrumMu_n->GetBinContent(ibin));
  }

  //---------------
  /*
  float vIntegral = hSpectrumMu->Integral();
  float integral = 0;
  float LowerLimitOfReachedMuon = 0; 
  for(int ibin=0;ibin<nbin;ibin++){
    integral += hSpectrumMu->GetBinContent(ibin+1);
    if( integral/vIntegral>(1-Ratio) ){LowerLimitOfReachedMuon=(ibin+1)*1.0e5/1.0e6; break;}
  }
  cout<<"Ratio = "<<Ratio<<endl;
  cout<<"Lower Limit of Muon Energy = "<<LowerLimitOfReachedMuon<<"MeV"<<endl;
  */
  /*------------------ Draw Histgram ----------------------*/
  gStyle->SetOptStat(000000);

  TCanvas *c1=new TCanvas("c1","c1",800,600);
  c1->SetGrid(1);
  c1->SetLogx(1);
  c1->SetLogy(1);

  hSpectrumMu->SetXTitle("Energy (MeV)");
  hSpectrumMu->SetYTitle("Flux (1/cm^{2}/s)");

  hSpectrumMu_p->SetLineStyle(7);
  hSpectrumMu_n->SetLineStyle(7);

  hSpectrumMu_p->SetLineColor(2);
  hSpectrumMu_n->SetLineColor(4);
  hSpectrumMu  ->SetLineColor(1);

  hSpectrumMu  ->Draw();
  hSpectrumMu_p->Draw("same");
  hSpectrumMu_n->Draw("same");

  TLegend *leg=new TLegend(0.3,0.25,0.5,0.5);
  leg->SetFillStyle(0);
  leg->AddEntry(hSpectrumMu_p,"#mu^{+}","l");
  leg->AddEntry(hSpectrumMu_n,"#mu^{-}","l");
  leg->AddEntry(hSpectrumMu  ,"#mu^{+} + #mu^{-}","l");
  leg->Draw();

  cout<<hSpectrumMu  ->GetMean()<<endl;

  /*
  TLine *l=new TLine(LowerLimitOfReachedMuon,0,LowerLimitOfReachedMuon,8.0e-6);
  l->SetLineColor(6);
  l->Draw();
  */

}
